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Abstract 

We investigate spherical macroions in the strong Coulomb coupling regime 
within the primitive model in salt-free environment. We first show that the 
ground state of an isolated colloid is naturally overcharged by simple elec- 
trostatic arguments illustrated by the Gillespie rule. We furthermore demon- 
strate that in the strong Coulomb coupling this mechanism leads to ionized 
states and thus to long range attractions between like-charged spheres. We 
use molecular dynamics simulations to study in detail the counterion distri- 
bution for one and two highly charged colloids for the ground state as well 
as for finite temperatures. We compare our results in terms of a simple ver- 
sion of a Wigner crystal theory and find excellent qualitative and quantitative 
agreement. 

PACS numbers: 82.70.Dd, 61.20.Qg, 41.20.-q 



Typeset using REVTeX 



* email: messina@mpip-mainz.mpg.de 
^email: holm@mpip-mainz.mpg.de 

* email : k. kremer @mpip-mainz . mpg.de 



1 



I. INTRODUCTION 



Charged colloidal suspensions are often encountered in the everyday life (technology 
biology, medicine ...) and have an important practical impact ffl. In numerous application- 
oriented situations, electrostatic repulsion among colloids (macroions) is desired in order to 
obtain a stabilized suspension. Consequently the understanding of the electrostatic interac- 
tion in such systems is motivated by practical as well as theoretical interests. There is recent 
experimental evidence that the effective interaction between two like-charged spherical col- 
loids (in the presence of neutralizing salts) can be attractive in the presence of one or two 
glass walls [|2|-§J. This is in contrast with the classical work of Derjaguin, Landau, Verwey 
and Overbeek (DLVO) based on a linearized Poisson Boltzmann theory |5|||, which foresees 
only repulsive effective Coulomb forces between two like-charged spheres even in confined 
geometry. There are some indications that this attraction might be explainable in terms of 
hydrodynamic effects induced by the walls [0]. 

Already in the bulk case there have been disputes for a long time about the existence of 
long range attractive forces, triggered mainly by the observation of voids in colloidal solutions 
H There is no clear experimental and theoretical picture, either, and there have been 



speculations that the experiments observed phase coexistence. Recent theoretical |pT^ [T4 
and simulation investigations have shown the existence of short range attraction. 



In two short communications |19|j20|] , we demonstrated by molecular dynamics (MD) 
simulations, how a mechanism involving overcharged and undercharged spherical macroions 
could lead to a strong long range attraction between charged spheres. In this paper we give 
a more detailed account and elaborate on the physical mechanism responsible for charge 
inversion (overcharge). Why and how does a charged particle strongly "bind" electrostat- 
ically at its surface so many counterions that its net charge changes sign? We further will 
discuss the necessary ingredients to explain this phenomenon in terms of a simple Wigner 
crystal theory. Using this Ansatz we show that it is possible for a pair of colloids which 
are sufficiently different in charge density to have an ionized ground state. Both, the one 
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and two colloid treated in terms of analytical predictions and verifications by 

simulation. Of special interest are the energy barriers necessary to cross from a neutral pair 
to an ionized pair state. We finally demonstrate by explicit simulations that the described 
features survive also at finite temperature case. 

The paper is organized as follows: in Sec. II a simple model based on the Gillespie rule is 
proposed to understand charge inversion. Section III contains details of our MD simulation 
model. Section IV is devoted for the study of a single highly charged colloid. In Sec. V 
we investigate the situation where two colloids are present. Finally, in Sec. VI we outline a 
summary of the results. 



II. UNDERSTANDING OVERCHARGING VIA THE GILLESPIE RULE 

Here we propose a simple model solely based on electrostatic energy considerations in 
order to understand the phenomenon of charge inversion for strongly coupled systems. Be- 
cause of the analogy between a spherical macroion surrounded by counterions and an atom [i. 
e. nucleus + electrons], it turns out be fruitful to use classical pictures of atomic physics in 
order to gain comprehension of certain phenomena occurring in mesoscopic colloidal systems 
1^ , 20] . To study the possibility of overcharging a single macroion, we recall the Gillespie 



rule also known as the valence-shell electron-pair repulsion (VSEPR) theory |2l|j2^] which 
is well known in chemistry to predict the molecular geometry in covalent compounds. Note 
that originally this model has nothing to do with overcharge. Applying simple electrostat- 
ics one can compute that the ground state structure of two, three, four and five electrons 
disposed on a hard sphere corresponds to simple geometrical situations like those depicted 
in Fig. |l|. The electrons try to maximize their mutual distances which leads, for example, 
in the case of 3 and 4 electrons to equilateral triangular and tetrahedral arrangements. 

Now, we can apply this concept to a spherical colloid of radius a, central charge Z m = 
+2e, where e is the elementary charge, and N c monovalent counterions. By referring to Fig. 
[I], the neutral system corresponds to the case where two counterions are present, and the 
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three other cases (three, four and five counterions) correspond to non-neutral overcharged 
states. 

The total electrostatic energy E(N C ) is merely made up of two terms: i) an attractive 
term E aU (N c ) due to the attraction between the counterions and the central charge and ii) 
a repulsive term E rep (N c ) due to the repulsion among the counterions. The final expression 
for the electrostatic energy as a function of the number of counter ions reads 

E(N C ) = E att (N c ) + E rep (N c ) = ksT— [-N c Z m + f{9)\ , (1) 

(X 

where = e 2 / '(4:7ieoe r k B T) is the Bjerrum length and f{6) is the repulsive energy part 
which is solely a function of the topology (relative angles between counterions, such as a and 
j3 appearing in Fig. [I], which also depend on N c ] of the ground state figure. For the specific 
cases reported in Fig. [1|, the calculation of E(N C ), with 2 < iV c < 5, is straightforward 
and the corresponding energy values are given in Fig. |I]. One deduces that the maximally 
obtainable overcharging is -2e (i. e. 100%) around the central charge. That is, the excess 
counterions gain more energy by assuming a topological favorable configuration than by 
escaping to infinity, the simple reason of overcharge. Note the arguments for overcharging are 
independent of the Bjerrum length and of the sphere radius, which enter only as prefactors 
in Eq. ©. 

To safely use this above outlined model one has just to ensure that the counterion size 
is small enough to avoid excluded volume effects, which in practice is always true. The 
important message is that, from an energy point of view, a colloid always tends to be 
overcharged. Obviously, for high central charge, the direct computation of the electrostatic 
energy by using the exact equation ([!]) becomes extremely complicated. Therefore we resort 
to simulations for highly charged spheres. 

III. SIMULATION MODEL 

The system under consideration contains two types of spherical charges: (i) one or two 
macroion(s) with a bare central charge Q = —Z m e (with Z m > 0) and (ii) small counterions 
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of diameter a with charge q = +Z c e (with Z c = 2) to neutralize the whole system. All these 
ions are confined in an impermeable cell and the macroion(s) is (are) held fixed. 

The molecular dynamics (MD) technique employed here is similar to this used in previous 
studies ||19| , |2"0|j . In order to simulate a canonical ensemble, the motion of the counterions 
is coupled to a heat bath acting through a weak stochastic force W(t). The equation of 
motion of counterion i reads 

m^ = -V i U-mry^ + W l (t), (2) 

where m is the counterion mass, U is the potential force having two contributions: the 
Coulomb interaction and the excluded volume interaction, and 7 is the friction coefficient. 
Friction and stochastic force are linked by the dissipation-fluctuation theorem < Wi(t) • 
Wj(t') >= Qm , yk B T8ij5(t — t'). For the ground state simulations the fluctuation force is set 
to zero. 

Excluded volume interactions are taken into account with a purely repulsive Lennard- 
Jones potential given by 

12 / \6" 
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+ e LJl for r - r < r cut , 

(3) 

0, for r - r > r cut , 



where r = for the counterion-counterion interaction, r = 7a for the macroion-counterion 
interaction, r cut (= 2 1 / 6 <r) is the cutoff radius. This leads to an effective macroion radius a 
(a = ro + cr = 8a) corresponding physically to the macroion-counterion distance of closest 
approach. Energy and length units in our simulations are defined as e^j =k B T (with 
To = 298 K) and a = 3.57 A respectively. In the following we will set ksTo = 1, so that all 
energies are measured in those units, suppressing thereby all factors of k B T in our equations. 

The pair electrostatic interaction between any pair ij, where % and j denote either a 
macroion or a counterion, reads 

U coul (r) = l B ——-^- , (4) 
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where Zi represents the valence of the ions (counterion or macroion). Being essentially 
interested in the strong Coulomb coupling regime we choose the relative permittivity e r = 16, 
corresponding to a Bjerrum length of 10er, for the remaining of this paper. To avoid image 
charges complications, the permittivity e r is supposed to be identical within whole the cell 
(including the macroion) as well as outside the cell. Typical simulation parameters are 
gathered in Table |. 



IV. ONE MACROION CASE 

In this section, we focus on counterion distribution exclusively governed by energy min- 
imization, i. e. T — OK. The single spherical macroion is fixed to the center of the large 
outer spherical simulation cell (i. e. both spheres are concentric) of radius R = 40<j. This 
leads to a colloid volume fraction f m = a 3 /R 3 = 8 x 10~ 3 . In such a case correlations are 
maximal, and all the counterions lie on the surface of the spherical macroion. To avoid being 
trapped in metastable states, we systematically heated and cooled (10 cycles) the system 



and only kept the lowest energy state then obtained |[23|| . It turns out that for this type 
of repulsive potential (between counterions) no rough energy landscape appears and thus, 
the MD method is efficient to find the ground state. First, we checked that this method 
reproduces well the ground state energies and structures of the simple situations depicted 
in Fig. |l| 



A. Counterion distribution 

To characterize the counterion layer structure, we compute the counterion correlation 
function g(r) on the surface of the sphere, defined as: 

c 2 g(r) = J25(r - ri)5(r - rj), (5) 

where c = N/4na 2 is the surface counterion concentration (N being the number of coun- 
terions), r corresponds to the arc length on the sphere. Note that at zero temperature all 
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equilibrium configurations are identical, thus only one is required to obtain g(r). The pair 
distribution g(r) is normalized as follows 

c 2irrg(r)dr= (N c + n-l), (6) 
Jo 

where N c = Z m /Z c is the number of counterions in the neutral state and n is the number of 
overcharging counterions. Because of the finite size and the topology of the sphere, g(r) has 
a cut-off at Tea (=25.1 a) and a zero value there. More precisely one cannot state that the 
uncorrelated case corresponds to g(r) = 1 for the present finite system. Therefore at "large" 
distance the correlation function differs from the one obtained with an infinite planar object. 
Furthermore the absolute value oi g(r) cannot be directly compared to the one obtained with 
an infinite plane. 

Correlation functions for the structural charge Z m = 180 and for two states of charge, 
neutral (n = 0) and overcharged (n = 8), can be inspected in Fig. [| One remarks that both 
structures are very similar and highly ordered. A snapshot of the ground state structure of 
the neutral state (n = 0) is depicted in Fig. [5| A visual inspection gives an almost perfect 
triangular crystalline structure (see Fig. [3p. A closer look at Fig. ^ reveals that the g(r) 
of the overcharged state, containing eight more counterions than the neutral one, shows its 
first peak at some shorter distance compared to the g(r) of the neutral state, as is expected 
for denser systems. 

It is also interesting to know how the counterion-layer structure looks like when the sys- 
tem is brought to room temperature T . At non zero temperature, correlation functions are 
computed by averaging <K r — r i)<K r — r j) over 1000 independent equilibrium configu- 
rations which are statistically uncorrelated. Results are depicted in Fig. f| for Z m = 180 
and / m = 8x 10~ 3 . As expected the long-range counterion positional order is neatly weaker 
at room temperature than in the ground state case. Meanwhile, the structure remains very 
correlated and highly short-range ordered and therefore it is referred as a strongly correlated 
liquid (SCL) [Q. In terms of Coulomb coupling parameter PUPS] T = Z^l B /a cc , where a cc 



is the average distance between counterions, we have r ~ 13 for Z m = 180. 



B. Energy analysis 



As demonstrated in Sec. II, the spatial correlations are fundamental to obtain overcharge. 
Indeed, if we apply the same procedure and smear Z counterions onto the surface of the 
colloid of radius a, we obtain for the energy 



E = I 



B 



1 Z Z m Z 



2 a 



(7) 



The minimum is reached for Z = Z m , hence no overcharging occurs. 

To generalize results of Sec. II to higher central charges we have considered three 
macroionic charge Z m of values 50, 90 and 180 corresponding to a surface charge density of 
one elementary charge per 180, 100 and 50 A 2 , respectively. For a given macroion, we always 
start by adding the exact number of counterions N c to have an electro-neutral system. Once 
equilibrium of this system is reached, we add the first overcharging counterion and let the 
new non-neutral system relax, and we repeat this operation a given number of times. The 
electrostatic energy is computed by summing up the pairwise interactions of Eq. (^) over 
all pairs. 

The electrostatic energy as a function of the number of overcharging counterions n is 
displayed Fig. ^. We note that the maximal (critical) acceptance of n (4, 6 and 8) increases 
with the macroionic charge Z m (50, 90 and 180 respectively). Furthermore for fixed n, the 
gain in energy is always increasing with Z m . Also, for a given macroionic charge, the gain 
in energy between two successive overcharged states is decreasing with n. 



The results of Sec. |IV A| showed that in the ground state the counterions were highly 



ordered. Rouzina and Bloomfield [25| first stressed the special importance of these crystalline 
arrays for interactions of multivalent ions with DNA strands, and later Shklovskii ( |T2|j24]| 
and references therein) showed that the Wigner crystal (WC) theory can be applied to 
determine the interactions in strongly correlated systems. In two recent short contributions 
T9| , |20[| we showed that the overcharging curves obtained by simulations of the ground state, 



like Fig. H, can be simply explained by assuming that the energy e per counterion on 



S 



the surface of a macroion depends linearly on the inverse distance between them, hence is 
proportional to y/~N for fixed macroion area, where N is the total number of counterions 



on the surface p!9|j20|j26|] . This can be justified by the WC theory. The idea is that the 
counterions form an ordered lattice on the surface of a homogeneously charged background 
of opposite charge, which is also called a One Component Plasma (OCP) Each ion 

interacts in first approximation only with the oppositely charged background of its Wigner- 
Seitz (WS) cell p4}| , which can be approximated by a disk of radius h, which possesses the 
same area as the WS cell. Because we can assume the area of the WS cell to be evenly 
distributed among the N counterions on the sphere's surface A = 4ira 2 we find 

*h* = ± = c~\ (8) 

The electrostatic interaction energy eW of one counterion with the background of its WS 
cell can then be determined by 



e (h) 



-l B Zl [ 2nrc-dr = -2^l B Z^, (9) 
J o r 



hence is proportional to yfc, which proves our initial assumption. It is convenient to define 
£ = IbZ% and aS h > = 2^Jtt fa 3.54. For fixed macroion area we can then rewrite Eq. (§) as 

e W(JV) = "-Tjv/iV. (10) 



If one computes this value for an infinite plane, where the counterions form an exact triangu- 
lar lattice, and takes into account all interactions, one obtains the same form as in Eq. (^j), 
but the prefactor gets replaced by the numerical value a wc = 1.96 ||29|| . Although the 



value is almost a factor of two smaller than the simple hole picture suggests, the functional 
dependence on the concentration is still the same. 

Not knowing the precise value of a we can still use the simple scaling behavior with 
c to set up an equation to quantify the energy gain AE\ by adding the first overcharging 
counterion to the colloid. To keep the OCP neutral we imagine adding a homogeneous 



surface charge density of opposite charge (— j^) to the colloid [28fl. This ensures that the 
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background still neutralizes the incoming overcharging counterion and we can apply Eq. 
(Pp. To cancel our surface charge addition we add another homogeneous surface charge 
density of opposite sign This surface charge does not interact with the now neutral 
OCP, but adds a self-energy term of magnitude |- , so that the total energy difference for 
the first overcharging counterion reads as 



AE 1 = (N c + l)e(N c + 1) - N c e(N c ) + —. 



(11) 



By using Eq. fllCf ) this can be rewritten as 



AEl = ^^U Nc + 1 )V2_ N 3/2 



f 2a 



(12) 



Completely analogously one derives for the energy gain AE n for n overcharging counterions 



al r 



[N c + nf' 2 - N 3 J 2 



in 2 
aT 



(13) 



Equation (|i~3| ) can be seen as an approximation of the exact general expression Eq. ([[]), 
where the topological term f(8) is handled by assuming a perfect planar crystalline structure 
through Eqs. (|TT]-|T^). Using Eq. (|T3"D, where we determined the unknown a from the 
simulation data for AE\ via Eq. (|l|) we obtain a curve that matches the simulation data 
almost perfectly, compare Fig. The second term in Equation ( |13D also shows why the 
overcharging curves of Fig. [5] are shaped parabolically upwards for larger values of n. 

Using the measured value of a we can simply determine the maximally obtainable number 
n ma x of overcharging counterions by finding the stationary point of Eq. ([13]) with respect to 
n: 

1/2 



Tlr, 



9a 2 3a 
+ 



1 



9a 2 



(14) 



32tt 4^ 

The value of n max depends only on the number of counterions iV c and a. For large iV c Eq. 

as the low temperature limit 
of a a neutral system in the presence of salt. What we have shown is that the overcharging in 



reduces to n max j^V^c which was derived in Ref. 
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this limit has a pure electrostatic origin, namely it originates from the topological favorable 
arrangement of the ions around a central charge. In the following we will investigate the 
behavior of a on the surface charge density and on the radius of the macroion. 

We have performed simulations for various surface charge densities by keeping A fixed 
and changing Z m = 2N C in the range 2 up to 180. Results can be found in Table [TJ and in 
Fig. p. We observe that a is already for values of N c as small as two, where one can use 
the Gillespie rule to calculate the energy exactly, close the planar value a^ wc \ and actually 
oscillates around this value. For Nc > 50, one reaches a plateau of a = 1.86 ± 0.05. 

This value is about 5% smaller then the one predicted by WC theory, and is presumably 
due to the finite curvature of the sphere. For large values of the radius a we expect a to 
reach the planar limit. To see the rate of convergence we varied^ a at a fixed concentration 
c. The results can be found in Table [III] and Fig. 0. For our smallest value of a = 6a 
we find a = 1.91. For small a, which is equivalent to a small number of N c , we observe 
again a slight oscillatory behavior of a, whereas for our two largest values a = 80a and 160o" 
we find up to numerical uncertainties the planar result a = a = 1.96. Again we stress 
that the numerical value of a enters only as a pref actor into the equations which govern the 
overcharging, it does not change the qualitative behavior. 

One could wonder if the results presented above are still valid when the bare central 
charge of the colloid is replaced by small discrete ions lying on the macroion surface? In 
fact it has been shown that the energy of the overcharged state (Fig. |5[) for random discrete 
colloidal charge distribution is more or less quantitatively affected f31|j3^] depending on 
the the valence of the counterions. More precisely it was shown that the overcharge still 
persists and has a similar (for monovalent counterions) or quasi-identical (for multivalent 
counterions) behavior to the one depicted in Fig. |5|, and this, even if ionic pairing occurs 



between the counterions and the discrete colloidal charges pT],|32|, that is even when no 



Note that this is the only part of the paper where a 7^ 8a. 
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counterion WC is formed. 



C. Macroion-counterion interaction profile 

In this part, we study the interaction potential profile at T = OK between a neutral 
effective macroion [bare macroion + neutralizing counterions] and one excess overcharging 
counterion at a distance r from the colloid center. The profile is obtained by displacing 
adiabatically the excess overcharging counterion from infinity towards the macroion. We 
investigated the case of Z m = 2,4,6,8,10,32,50,90,128,180, and 288. All curves can be 
nicely fitted with an exponential fit of the form 

E^r) = AE ie ~ T{r - a \ (15) 

where AEi is the measured value for the first overcharging counterion, and r is the only fit 
parameter (see Table [II]). Results for the two values Z m = 50 and 180 are depicted in Fig. 
|8]. If one plots all our results for r versus y/N~ c we observe a linear dependence for a wide 
range of values for iV c , 



T = m^/N c , (16) 

with ma ~ 0.1, as can be inspected in Fig. || 

This behavior can again be explained using a "WC hole" picture in the limiting situation 
where x := r — a is small (i. e. the displaced counterion is close to the macroion surface). 
To this end we consider the classical electrostatic interaction V<n s k{ x ) between a uniformly 
charged disk (the WC hole - supposed planar) and a point ion (the displaced counterion) 
located on the axis of the disk at a distance x from its surface, which is given by 



V disk {x) = -2n£c(Vh? + x^-x). (17) 



As in Eq. (^), h = (ire) x l 2 is the hole radius. For small distance x, we expand Eq. (|T 



V disk {x) = e ih) 



r 
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(18) 



where the surface term Vdi s k{ x = 0) = is given by Eq. @. By expanding the exponential 
in Eq. fll5|) to 2 nd order for small rx we obtain 



E 1 {x)=AE 1 



1-tx + —x 2 + 0{t 3 x 3 ) 



(19) 



A comparison between Eq. (|T9| ) and Eq. (|T8| ) shows that to this order we can identify 



1 

h 



7TC 



2a 



0.06a/ N r . 



(20) 



Comparing this to Eq. (|16D we note that this simple illustration gives us already the correct 
scaling as well as the prefactor up to 30%. We neglected here the effect that the surface 
concentration changes when the ion is close to the macroion as well as the curvature of the 



macroion. 



V. TWO MACROIONS CASE 

In this section we consider two fixed charged spheres of bare charge Qa and Q b separated 
by a center-center separation R and surrounded by their neutralizing counterions. All these 
ions making up the system are immersed in a cubic box of length L = 80a, and the two 
macroions are held fixed and disposed symmetrically along the axis passing by the two centers 
of opposite faces. This leads to a colloid volume fraction f m = 2 • |7r(a/L) 3 « 8.4 x 10~ 3 . 
For finite colloidal volume fraction f m and temperature, we know from the study carried 
out above that in the strong Coulomb coupling regime all counterions are located in a 
spherical "monolayer" in contact with the macroion. Here, we investigate the mechanism of 
strong long range attraction stemming from monopole contributions: that is one colloid is 
overcharged and the other one undercharged. 



A. Like charged colloids 
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1. Observation of metastable ionized states 

In the present charge symmetrical situation we have Qa = Qb = —Z m e. This system 
is brought at room temperature T . Initially the counterions are randomly generated inside 
the box. Figure [H] shows two macroions of bare charge Z m = 180 surrounded by their 
quasi-two-dimensional counterions layer. The striking peculiarity in this configuration is 
that it corresponds to an overcharged and an undercharged sphere. There is one counterion 
more on the left sphere and one less on the right sphere compared to the bare colloid charge. 
Such a configuration is referred as ionized state. In a total of 10 typical runs, we observe 
this phenomenon 5 times. We have also carefully checked against a situation with periodic 
boundary conditions, yielding identical results. However it is clear that such a state is 
"metastable" because it is not the lowest energy state. Indeed, in this symmetrical situation 
the ground state should also be symmetrical so that both colloid should be exactly charge- 
compensated. Such arguments remain valid even at non-zero temperature as long as the 
system is strongly energy dominated, which is presently the case. Nevertheless the ionized 
states observed here seem to have a long life time since even after 10 8 MD time steps this 
state survives. In fact we could not observe within the actual computation power the recover 
of the stable neutral state. To understand this phenomenon we are going to estimate the 
energy barrier involved in such a process. 

2. Energy barrier and metastability 

To estimate the energy barrier, electrostatic energy profiles at zero temperature were 
computed, where we move one counterion from the overcharged macroion to the under- 
charged, restoring the neutral state [see drawing depicted in Fig. [II|(a)1 . We have checked 
that the path leading to the lowest barrier of such a process corresponds to the line joining 
the two macroions centers. The simulation data are sketched in Figs. |Tl~|(a-b) and were fitted 
using a similar technique to the single macroion-counterion interaction profile given by Eq. 
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(|15l), which will be explicitly treated later. The resulting simulated energy barrier AEb ar 
is obtained by taking the difference between the highest energy value of the profile and the 
ionized state energy (start configuration). Values of AEt, ar can be found in Table [IV] for the 
small macroion separation case R/a = 2.4. One clearly observes a barrier, which increases 
quasi linearly with the charge Z m for the small colloids separation R/a = 2.4 [cf. Fig. |TT|(a) 
and Table |V|]. The ground state corresponds as expected to the neutral state. Note that 
the ionized state and the neutral state are separated by only a small energy amount (less 
than 2.5), the difference being approximately of the order of the monopole contribution 
E = Z#(4/8 — 4/11) « 1.36. The physical origin of this barrier can be understood from 
the single macroion case where we showed that a counterion gains high correlational energy 
near the surface. This gain is roughly equal for both macroion surfaces and decreases rapidly 
with increasing distance from the surfaces, leading to the energy barrier with its maximum 
near the midpoint. For the single macroion case we showed that the correlational energy 
gain scales with y/ Z m , whereas here we observe a linear behavior of the barrier height with 
Z m . We attribute this effect to additional ionic correlations since both macro ions are close 
enough for their surface ions to interact strongly. For large separations (here R/a = 4.25) 
we find again that the barrier height increases with y/Z m) as expected [see Fig. |TT](b) and 
Table [V]]- Furthermore the energy barrier height naturally increases with larger colloidal 
separation. The Z m dependence of the barrier also shows that at room temperature such 
ionized states only can occur for large Z m . In our case only for Z m = 180, the ionized state 
was stable for all accessible computation times. Unfortunately, it is not possible to get a 
satisfactory accuracy of the energy jumps at non-zero temperatures. Nevertheless, since we 
are interested in the strong Coulomb coupling regime, which is energy dominated, the zero 
temperature analysis is sufficient to capture the essential physics. 

Simulation results presented in Fig. [11| can be again theoretically well described using the 
previously exponential profiles obtained for the macroion-displaced counterion in Sec. [IV C 



for a single colloid. For the two macroions case, the general expression for the electrostatic 
interaction Ebar{r,R) of the present process can be approximated as 
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E bar (r, R) = AE{ exp [-r(r - a)] + AE{ exp [-r(R - r - a)] - , (21) 

where AS* is the "effective" correlational energy gained by the first OC at one macroion 
surface assumed identical for both colloids. The last term in Eq. ( pT| ) corresponds to the ad- 
ditional monopole attractive contribution of the displaced counterion with the undercharged 
colloid. Fitting parameters (AEl an d t) for R/a = 2.4 and R/a = 4.25 can be found in Ta- 
bles [IV] and |VJ respectively. Same values of r were used here as those of the single macroion 
case (see Fig. |S] and Table [TJ) . However for the small colloidal separation (R/a = 2.4), 
due to the extra inter-colloidal surface counterions correlations, we used a slightly larger 
(absolute) value for AEl compared to the one (AEi) of an isolated colloid (compare Ta- 
bles |V| and |V|) . This is compatible with the idea that between the two colloids (especially 
when both spheres come at contact), we have the formation of a "super-layer" which is 
more dense, thus leading to a smaller hole radius and a higher energy gain. An analysis of 
the counterions structure of the two macroions reveals that both WC counterion layers are 
interlocked, that is the projection along the axis passing through the colloid centers gives a 
superlattice structure (see Fig. [12]). 

For large colloidal separation (R/a = 4.25), the WC structure on one of the colloids is 
unperturbed by the presence of the other, hence we can take AEl = AEi, and our simulation 
data can nicely be fitted by the parameters inferred from the single colloid system. 



3. Effective forces 

Results concerning the effective forces at zero temperature between the two macroions 
are now investigated which expression is given by 

Feff(R) = F mm (R) + F LJ + F mc , (22) 

where F mm (R) is the direct Coulomb force between macroions, Flj is the excluded volume 
force between a given macroion and its surrounding counterions and F mc is the Coulomb 
force between a given macroion and all the counterions. Because of symmetry, we focus 
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on one macroion. To understand the extra-attraction effect of these ionized-like states, we 
consider three cases: (i) F ion = F e ff in the ionized state with a charge asymmetry of ± 1 
counterion (ii) F neut = F e ff in the neutral case (iii) F mono = F e ff simply from the effective 
monopole contribution. Our results are displayed in Fig. [13] for Z m = 180 , where the 
ionized state was also observed at room temperature. The non-compensated case leads to 
a very important extra attraction. This becomes drastic for the charge asymmetry of ± 
2 counterions at short separation R/a = 2.4 leading to a reduced effective attractive force 
FIb = —10.7, a situation which was also observed in our simulation at room temperature. In 
contrast to previous studies fl~5| , |l6|l , these attractions are long range. For a sufficiently large 
macroion separation (from 3.5a), corresponding here roughly to a macroion surface-surface 
separation of one colloid diameter, the effective force approaches in good approximation the 



monopole contribution (see Fig. [T3|). 



B. Asymmetrically charged colloids 

In this section we investigate the case where the two colloids have different charge den- 
sities. We will keep the colloidal radii a fixed, but vary the bare colloidal charges. The 
charge on sphere A is fixed at Za = 180, and sphere B carries variable charges with Zb 
(where Zb < Za) ranging from 30 up to 150. Global electroneutrality is ensured by adding 
Na + Nb divalent counterions, with Na = Za/Z c , and Nb = Zb/Z c . In this way we vary 
the bare counterion concentrations = JNi/Aira 2 , where i stands for A or B. 



1. Ground state analysis 

We start out again with studying the ground state of such a system. The electrostatic 
energy of the system is investigated for different uncompensated bare charge cases (ionized 
states) by simply summing up Eq. (|]) over all Coulomb pairs. We define the degree of 
ionization (DI) as the number of counterions overcharging colloid A (or, equivalently, un- 
dercharging colloid B). The system is prepared in various DI and we measure the respective 
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energies. These states are separated by kinetic energy barriers, as was demonstrated above. 
We consider three typical macroionic charges Zb (30, 90 and 150) and separations R/a (2.4, 



3.0 and 4.25). The main results of the present section are given in Fig. [TJ]. For the largest 
separation R/a = 4.25 and largest charge Zb = 150 [see Fig. 0(a)], one notices that the 
ground state corresponds to the classical compensated bare charge situation [referred as the 
neutral state (DI=0)}. Moreover the energy increases stronger than linear with the degree 
of ionization. If one diminishes the bare charge Zb to 90 and 30, the ground state is actually 
the ionized state for a DI of 1 and 3, respectively. The ionized ground state is about 8 
and 36 , respectively, lower in energy compared to the neutral state. This shows that even 
for a relative large colloid separation, stable ionized states should exist for sufficient low 
temperatures and that their stability is a function of their charge asymmetry. 

For a shorter separation R/a = 3.0, ionized ground states are found [see Fig. 0(b)] for 
the same charges Zb as previously. Nevertheless, in the ground state the DI is now increased 
and it corresponds to 2 and 4 for Zb = 90 and 30 respectively. The gain in energy is also 
significantly enhanced. For the shortest separation under consideration R/a = 2.4 [see Fig. 
|H|(c)], the ground state corresponds for all investigated values of Z B to the ionized state, 
even for Zb = 150. We conclude that decreasing the macroion separation R enhances the 
degree of ionization and the stability of the ionized state. 

To understand this ionization phenomenon, it is sufficient to refer to an isolated macroion 
surrounded by its neutralizing counterions. We have investigated the energies involved 
in the ionization (taking out counterions). The complementary process of overcharging 
(adding counterions) has already been investigated (see Fig. |^). A derivation of the formula 
describing the ionization energy AE lon proceeds completely analogously to the one carried 
out for the overcharging Eq. flT3| ) and gives for the n th degree of ionization 



Bn _ i) ^2 



{NB-n)^-NT]+[\, (23) 



'A 

where a A ' B are the values of a belonging to colloid A and B, respectively. In Fig. we 



compare the predictions of Eqs. fll3| , £3|) to our simulation data, which shows excellent 



agreement. Our numerical data for AE[ on for Nb = 15, 45, and 75, the value of AE± C for 
Na = 90 (overcharging process), as well as the corresponding values for a, which have been 
used for Fig. [15] can be found in Table [TT]. 



With the help of Eqs. (|T3| , |23|) , one can try to predict the curves of Fig. [H] for finite 
center-center separation R. Using for colloid A and B the measured values ot A and a B , we 
obtain for the electrostatic energy difference at finite center-center separation R 
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The quality of the theoretical curves can be inspected in Fig. 0. The prediction is 
is very good for large separations, but the discrepancies become larger for smaller separa- 
tions, and one observes that the actual simulated energies are lower. Improvements could 
be achieved by including polarization effects along the ideas leading to Eq. ((|21~1), by ad- 
justing, for example, a A and a B . More important, the physical interpretation of Eq. (|24|) 
is straightforward. The left two terms represent the difference in correlation energy, and 
last term on the right the monopole penalty due to the ionization and overcharging process. 
This means that the correlational energy gained by overcharging the highly charged colloid 
A must overcome the loss of correlation energy as well as the monopole contribution {two 
penalties) involved in the ionization of colloid B. With the help of Eq. (|24] ) we can establish 
a simple criterion (more specifically a sufficient condition), valid for large macroionic sepa- 
rations, for the charge asymmetry \/Na — V^b to produce an ionized ground state of two 
unlike charged colloids with the same size: 



'Nb) > 



3a 



(25) 



Referring to Fig. (|15|) this criterion is met when the overcharge curve (changed sign) is 
higher than the ionization curve. 
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If one uses the parameters of the present study one finds the requirement Nb < 66 to get 



a stable ionized state. This is consistent with our findings where we show in Fig. [1^ that for 
Nb = 75, and R/a = 4.25, no ionized ground state exists whereas for Nb = 60 we observed 
one even for infinite separation (not reported here). The criterion Eq. (p5| ) is merely a 
sufficient condition, since we showed in Fig. [TJ] that when the colloids are close enough 
this ionized state can appear even for smaller macroion charge asymmetry due to enhanced 
inter-colloidal correlations. At this stage, we would like to stress again, that the appearance 
of a stable ionized ground state is due merely to correlation. An analogous consideration 
with smeared out counterion distributions along the lines of Eq. ([^D will again always lead 
to two colloids exactly neutralized by their counterions p3| |. Our energetical arguments are 



quite different from the situation encountered at finite temperatures, because in this case 
even a Poisson-Boltzmann description would lead to an asymmetric counterion distribution. 
However, in the latter case this happens due to pure entropic reasons, namely in the limit 
of high temperatures, the counterions want to be evenly distributed in space, leading to an 
effective charge asymmetry. 

At this stage, on looking at the results presented above, it appears natural and straight- 
forward to establish an analogy with the concept of ionic bonding. It is well known in 
chemistry that the electro-negativity concept provides a simple yet powerful way to predict 
the nature of the chemical bonding |34j] . If one refers to the original definition of the electro- 



negativity given by Pauling p4| : "the power of an atom in a molecule to attract electrons 
to itself" , the role of the bare charge asymmetry becomes obvious. Indeed, it has an equiv- 
alent role at the mesoscopic scale as the electron affinity at the microscopic scale. Another 
interesting analogy is the influence of the colloidal separation on the stability of the ionized 
state. Like in diatomic molecules, the ionized state will be (very) stable only for sufficiently 
short colloid separations. Nevertheless, one should not push this analogy too far. One point 
where it breaks down concerns the existence of an ionized ground state in colloidal system 
for large colloid separation, providing that the difference in the counterion concentration 
on the surface is large enough. In an atomistic system this is impossible since even for the 
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most favorable thermodynamical case, namely CsCl, there is a cost in energy to transfer an 
electron from a cesium atom to a chlorine atom. Indeed, the smallest existing ionization 
energy (for Cs, 376 kJ mol -1 ) is greater in magnitude than the largest existing electron 
affinity (for Cs, 349 kJ mol -1 ). In other terms, for atoms separated by large distances in 
the gas phase, electron transfer to form ions is always energetically unfavorable. 



2. Finite temperature analysis 

As a last result, aimed at experimental verification, we show that an ionized state can 
also exist spontaneously at room temperature T . Figure [16| shows the time evolution of the 
electrostatic energy of a system Za = 180 with Zb = 30, R/a = 2.4 and f m = 7 ■ 10 -3 , 
where the starting configuration is the neutral state (DI = 0). One clearly observes two 
jumps in energy, AEi = —19.5 and AE 2 = — 17 A , which corresponds each to a counterion 
transfer from colloid B to colloid A. These values are consistent with the ones obtained for 
the ground state, which are— 20.1 and —16.3 respectively. Note that this ionized state (DI 
= 2) is more stable than the neutral but is expected to be metastable, since it was shown 
previously that the most stable ground state corresponds to DI = 5. The other stable 
ionized states for higher DI are not accessible with reasonable computer time because of 
the high energy barrier made up of the correlational term and the monopole term which 



increases with DI . In Fig. [17| we display a typical snapshot of the ionized state (DI = 2) 
of this system at room temperature. 

Obviously, these results are not expected by a DLVO theory even in the asymmetric case 
(see e. g. ||35|| ). Previous simulations of asymmetric (charge and size) spherical macroions 
were also far away to predict such a phenomenon since the Coulomb coupling was weak 



(water, monovalent counterions) . 
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VI. CONCLUDING REMARKS 



In summary, we have shown that the ground state of a charged sphere in the presence 
of excess counterions is always overcharged. A sufficiently charged colloid can in principle 
be highly overcharged due to counterion mediated correlation effects, and this phenomenon 
is quantitatively well described by a simple version of Wigner crystal theory In the strong 
Coulomb coupling regime, the energy gain of a single excess ion close to a counterion layer 
can be of the order of many tens of /c_bT . Furthermore we demonstrated that the electrostatic 
interaction between a counterion and a macroion effectively neutralized by its counterions 
decays exponentially on a length scale which is equal to the Wigner crystal hole radius. 

We further found that for two like-charged macroions (symmetric case), an initially 
randomly placed counterion cloud of their neutralizing divalent counterions may not be 
equally distributed after relaxation, leading to two macroions of opposite net charges. This is 
due to the short range WC attraction which leads to this energetically favorable overcharged 
state. The resulting configuration is metastable, however separated by an energy barrier of 
several UbTq when the bare charge is sufficiently large, and can thus survive for long times. 
Such configuration possess a natural strong long range attraction. 

In return, if the symmetry in the counterion concentration on the colloidal surface is 
sufficiently broken, the ionized state can be stable. The ground state of such a system is 
mainly governed by two important parameters, namely the asymmetry in the counterion 
concentration determined by y/cX — \fc~Bi and the colloid separation R. If the counterion 
concentration difference is high enough, the ground state corresponds to an ionized state, 
whatever the macroions separation R is. However, the degree of ionization depends on R. 
Besides, for large R, we have established a criterion, allowing to predict when a stable ionized 
configuration can be expected. The counterion concentration difference plays an analogous 
role to the electron affinity between two atoms forming a molecule with ionic bonding. We 
demonstrated that the results presented here for the ground state can lead to a stable ionic 
state even at room temperature providing that the Coulomb coupling and/or the counterion 
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concentration asymmetry is sufficiently large. This is also a possible mechanism which could 
lead to strong long range attractions, even in bulk. Future work will treat the case where 
salt ions are present. 
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TABLES 



TABLE I. Simulation parameters with some fixed values. 



parameters 

a = 3.57 A Lennard Jones length units 

To = 298K room temperature 

€lj = ksTo Lennard Jones energy units 

Z m macroion valence 

Z c = 2 counterion valence 

Ib = lOd Bjerrum length 

f m macroion volume fraction 
a = 8a macroion-counterion distance of closest approach 
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TABLE II. Measured values for an isolated macroion, with fixed radius a, of the energy gain for 
the first overcharging counterion AE® C for various macroion bare charge Z m = 2N C . The value of 
a can be compared to the prediction of WC theory for an infinite plane, which gives 1.96, compare 
text. We also record the values of the fitting parameter r of Eq. (|l~5| ) for selected N c corresponding 
to those of Fig. (||). The symbol " stands for the ionization process discussed in Sec. VB 1 . 





N c 




a 


TO 


2 


1 


-2.5 


1.94 


0.12 


4 


2 


-3.8 


1.89 


0.18 


6 


3 


-5.3 


1.97 


0.19 


8 


4 


-6.1 


1.92 


0.24 


10 


5 


-7.5 


2.02 


0.24 


20 


10 


-10.7 


1.93 


- 


30« 


15 


+17.9 


1.91 


- 


32 


16 


- 


- 


0.41 


50 


25 


-18.0 


1.92 


0.51 


90 


45 


-24.4 


1.88 


0.68 


90W 


45 


+29.2 


1.89 




128 


64 






0.79 


150® 


75 


+37.4 


1.91 




180 


90 


-35.3 


1.88 


0.93 


288 


144 






1.19 


360 


180 


-50.0 


1.86 
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TABLE III. Measured values of the energy gain AE? c and fixed counterion concentration c, 



varyinj 


I this time the macroion radius a and the number of counterions N c . 




a/a 


N c 


AE 1 /k B T Q 


a 


6 


9 


-13.3 


1.91 


8 


16 


-14.4 


1.97 


10 


25 


-14.5 


1.93 


12 


36 


-14.7 


1.92 


14 


49 


-15.1 


1.94 


16 


64 


i r i 

-15.1 


1.92 


20 


100 


-15.3 


1.92 


40 


400 


-15.9 


1.94 


80 


1600 


-16.4 


1.97 


160 


6400 


-16.5 


1.96 


TABLE IV. Measured value of the energy 


barrier and fit parameters of the electrostatic inter- 


action 


process involved in Fig. |ll](a) for R/a ■ 


= 2.4 and for for different macroion bare charj 


;es. 




AE bar /k B T 


AE$/k B T 


TO 


50 


4.9 


-20.4 


0.51 


90 


9.6 


-27.5 


0.68 


180 


20.4 


-39.4 


0.92 
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TABLE V. Measured value of the energy barrier and fit parameters of the electrostatic inter- 
action process involved in Fig. |ll](b) for R/a = 4.25 and for different macroion bare charges. 

Z m AE bar /k B T AEl/keTo to 

50 16.8 -18.4 0.51 

90 23.3 -24.4 0.68 

180 33.8 -35.3 0.92 
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FIGURES 

FIG. 1. Ground state configurations for two, three, four and five counterions. The correspond- 
ing geometrical figures show the typical angles. The electrostatic energy (in units of ksTls/a) is 
given for a central charge of +2e. 

FIG. 2. Ground state surface counterion correlation functions for Z m = 180 and two states of 
charge [neutral (n = 0) and overcharged (n = 8)]. 

FIG. 3. Snapshot of the ground state structure of the neutral state (n = 0) with a macroion 
charge Z m = 180 [see Fig. [2] for the corresponding g(r)]. 

FIG. 4. Surface counterion correlation functions at room temperature Tq for two states of charge 
[neutral (n=0) and overcharged (n=8)] with Z m = 180 and f m = 6.6 x 10~ 3 . 

FIG. 5. Electrostatic energy (in units of HbTq) for ground state configurations of a single 
charged macroion of as a function of the number of overcharging counterions n for three different 
bare charges Z m . The neutral case was chosen as the potential energy origin, and the curves were 
produced using the theory of Eq. (|l3|), compare text. 



FIG. 6. Wigner crystal parameter a as a function of the number of counterions iV c for fixed 
colloid radius a. 

FIG. 7. Wigner crystal parameter a as a function of the colloid radius a for a fixed surface 
counterion concentration c. 

FIG. 8. Electrostatic interaction energy (in units of &bTo) of a divalent counterion with a 
neutral effective colloid [bare particle + surrounded counterions] as function of distance r/a from 
the center of a macroion for two different macroion bare charges Z m . The energy is set to zero at 



distance infinity. Solid lines correspond to exponential fits [see Eq. (15)]. 



FIG. 9. Exponential fit parameter r as a function of the square root of the number of counte- 
rions y/Nc. The dashed line corresponds to a linear fit in \fN~ c - 
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FIG. 10. Snapshot of a "pseudo-equilibrium" configuration at room temperature To where the 
counterion-layers do not exactly compensate the macroions charge. Here the deficiency charge is 
±1 counterion (or ±2e as indicated above the macroions) and R/a = 3.6. 

FIG. 11. Total electrostatic energy (in units of &,bTo) of the system, for zero temperature 
configurations, of two macroions at a center-center separation of (a) R/a = 2.4 (b) R/a = 4.25 
as a function of one displaced counterion distance from the left macroion for three typical values 
Z m . The exact neutral state was chosen as the potential energy origin. The schematic drawing 
indicates the path (dotted line) of the moved counterion. The ending arrows of the arc indicate 
the start position (left sphere) and final position (right sphere) of the moved counterion. Dashed 
lines correspond to the fit using Eq. (|2l|) of which parameters can be found in Tables [V and [V|. 



FIG. 12. Projection of the counterion positions, located on both inner (face to face) hemi- 
spheres, along the symmetrical axis passing through the macroion centers. Open (filled) circles 
are counterions belonging to macroion A (B). One clearly sees the interlocking of the two ordered 
structures yielding locally to a superlattice. 

FIG. 13. Reduced effective force between the two spherical macroions at zero temperature for 
Z m = 180 as a function of distance from the center. The different forces are explained in the text. 
The lines are a guide to the eye. 

FIG. 14. Total electrostatic energy as a function of the degree of ionization for zero temperature 
configurations of two colloids [A and B), for three typical charges Zb (30, 90 and 150) for macroion 
B and for three given distance separations: (a) R/a = 4.25, (b) R/a = 3.0 and (c) R/a = 2.4. 
Dashed lines were obtained using Eq. (|2j) . 
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FIG. 15. Total electrostatic energy as a function of the degree of ionization for zero temperature 
configurations of an isolated colloid. The three upper curves correspond to the ionization energy 
for the three typical charges Zb (30, 90 and 150). The lower curve corresponds to the energy 
gained (changed sign for commodity) by overcharging (Za = 180). Dashed lines were obtained 



using Eqs. (O, 23) with the measured values for a from Table O. 



FIG. 16. Relaxation, at room temperature To = 298-ftT, of an initial unstable neutral state 
towards ionized state. Plotted is the total electrostatic energy versus time (LJ units), for Zb = 30 
and R/a = 2.4. Dashed lines lines represent the mean energy for each DI state. Each jump in 
energy corresponds to a counterion transfer from the macroion B to macroion A leading to an 
ionized state (DI = 2) which is lower in energy than the neutral one. The two energy jumps 
AEi/ksTQ = —19.5 and AE^/^bTo = —17.4 are in very good agreement with those of Fig. |i~4|(c) 
(-20.1 and -16.3). 

FIG. 17. Snapshot of the ionized state (DI = 2) obtained in the relaxation process depicted in 
Fig. [l(| with the net charges +4e and — 4e as indicated. 
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